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Abstract 

Using the dynamical matrix of a crystal obtained from ab initio calculations, we have evaluated 
for the first time the strength of the dynamic flexoelectric effect and found it comparable to that of 
the static bulk flexoelectric effect, in agreement with earlier order-of-magnitude estimates. We also 
proposed a method of evaluation of these effects directly from the simulated phonon spectra. This 
method can also be applied to the analysis of the experimental phonon spectra, being currently 
the only one enabling experimental characterization of the static bulk flexoelectric effect. 
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The flexoelectric effect, which is a polarization response to a strain gradient, has been at¬ 
tracting appreciable attention of theorists and experimentalists. Being a higher-order effect 
with respect to piezoelectricity, it becomes appreciable at the nanoscale, where large strain 
gradients are expected. Recent experimental studies attest to many flexoelectricity-driven 
phenomena. For example, a strain gradient can work as an electric held via hexoelectric 
coupling: it can induce poling, switching, and rotation of polarization nm. It can create 
a voltage offset of hysteresis loops and smear the dielectric anomaly at ferroelectric phase 
transitions [3]. The hexoelectric ehect actually represents 4 related phenomena: static and 
dynamic bulk hexoelectric ehects, surface hexoelectric ehect and surface piezoelectricity [B]. 
Among the four ehects only the static bulk hexoelectric ehect can be viewed as a high-order 
analog of piezoelectricity. It is the only ehect that has been quantitatively assessed by the¬ 
orists [MU- In contrast, dynamic hexoelectricity has never been quantitatively assessed 
either theoretically or experimentally. Though, a general theory of the ehect has being 
ohered [12]. Concerning the size of the ehect, it was only demonstrated that the dynamic 
ehect is expected to be comparable to the static bulk hexoelectricity [B] . Thus, one raises the 
reasonable question of the real size of dynamic hexoelectricity. It is an important question 
since nowadays papers on electromechanical modelling involving hexoelectricity always ne¬ 
glect the dynamic hexoelectric ehect. Another issue missing from the current state-of-the-art 
of hexoelectricity is a method of experimental evaluation of this ehect. 

In this Letter we present the results of quantitative modeling of the dynamic hexoelectric 
ehect in cubic SrTiOs (STO), the classical perovskite material with the most studied hexo¬ 
electric properties. We will demonstrate that the size of the dynamic hexoelectric ehect is 
indeed comparable to that of the static bulk hexoelectricity, identifying a situation where 
the dynamic ehect is a few times larger than the static. We will also oher a method for the 
evaluation of the dynamic hexoelectric ehect from simulated phonon spectra; the method 
that can also be used for assessment of this ehect using experimental phonon data. 

In terms of electromechanical constitutive equations |B| 

Ei = X-^Pj - + 7 ./;, ( 1 ) 


ljUJbl 

duki d^Pk 

pUi = Cijki^ -h - MjiPj^ 


dxj ' dxidxi 

the static and dynamic bulk hexoelectric ehects are introduced via the flexocoupling tensor 
fkiij and tensor Mjj, which will be termed as flexodymanic tensor^ respectively. Equations 
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Q and (|^ are written in the Cartesian reference frame Xi, summation over repeating suffixes 
is accepted. We use the following notation: Pj - polarization, Ej - macroscopic electric held, 
Uj - elastic displacements, Uij = ^{dUi/dxj + dUj/dxi) - strain tensor in the approximation 
of small strains. 

When one is interested in the macroscopic (static or dynamic) response, one can neglect 
the last two terms in Eqs. ([^ and ([^, which correspond to higher dispersions. Then, for 
the case of the static strain gradient {Uj = 0), dehning the hexoelectric response under the 
condition of a vanishing macroscopic held [B] , Eq. ([^ yields 

Pk = (3) 

OXi 

where the hexoelectric tensor reads as 


f^ijkl X,kxfijxl- (4) 

Equations (|^ and (|^ describe the static bulk hexoelectric ehect. 

In the dynamic case, eliminating Uj between ([^ and (|^ one hnds an alternative expres¬ 
sion for the hexoelectric tensor 

fJ^ijkl = Xkxfijxl- ( 5 ) 


fijxl = fijxl - P MxmCmlij- (6) 

Equation ([^ introduces the total flexocoupling tensor fU^i, where the term containing 
p~^MxmCmiij coutrols the contribution of the dynamic hexoelectricity to the total bulk hex¬ 
oelectric response. As is clear from Eq. ([^ the origin of the dynamic ehect is an accelerated 
motion of the medium. Since in an elastic wave Uj oc the corresponding polarization 
response is proportional to the strain gradient. 

In this work we are interested in the evaluation of the dynamic bulk hexoelectric ehect and 
its comparison with the static ehect for an STO crystal, in the Born-charge approximation. 
To do this, the tensors M^m and fij^i were obtained using the expansion of the dynamical 
matrix Aipj/pi{q) calculated with the long-range contribution of the macroscopic electric 
held being excluded na. Here, the suffixes i and i' enumerate the Cartesian components of 
the atomic displacements while p and p' - atoms in an elementary unit cell of the crystal; 
q denotes the wave-vector. Specihcally, one needs the coefficients of the following long- 
wavelength expansion of the dynamical matrix 

Aipj'p'i^ = + Afjl,^,qj + -AfJj,^p,qjqi + ..., 


(7) 
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and the matrix Fip^j/p/, that is the inverse of the singnlar matrix dehned in a special 

way [13]. The resnlts of Tagantsev’s approach [T2|, written for a cnbic crystal of perovskite 
structnre in the Born charge approximation, can be summarised as follows 


Mij 


X 


-1 


Qix,p^xp,jp' ij^p' 




fm = - K%), 


-1 


Ajjkl _ ( A2)kl Vp" ^(2)kl . 

■^'tp 9 / A wXp'.iv" / j .ia")' 


( 8 ) 

(9) 

( 10 ) 


p/f Q'Q'' 

where all summations are done over all s atoms of the unite cell. Here v is the volume of 
the unit cell; Qij^p and rup are the matrix of the Born charges and the mass of p-th atom in 
the unit cell; y is the only non-zero component (in a cubic crystal) of the tensor of dielectric 
susceptibility [T^ 

Xij ~Qix,p^xp,yp'Qjy,p'• (1^) 

Thus to hnd the tensors Xxmi and fijxi, only the matrices and Af^^!^, (available 

once the dynamical matrix jy (g) is calculated), the matrix of the Born charges Qix,p and 
the volume of the unit cell are needed. 

We found Aip^iipi{q) and Qix,p using hrst principles calculations. Since, within zero Kelvin 
DFT calculations, the cubic structure of STO is unstable (negative squares of phonon fre¬ 
quencies for the soft-mode phonons) we performed our calculations under a hydrostatic 
pressure of 78 kBar to stabilize it, corresponding to approximately 1% elastic strain. Since 
the tensors Mxm and fijxi are not critically small, in view of the Landau theory arguments 
and atomic estimates, such strain is expected to induce some 1% modihcations of these 
tensors. On the same lines, the result of zero-Kevin calculations should give a reasonable 
approximation of the hnite-temperature values of these tensors. Such an approach is justi- 
hed by the experience of zero-Kevin calculations of elastic constants. This way, we believe 
that the estimates obtained here are valid for realistic hnite temperature and pressure. 

The evaluation of the contribution of dynamic flexoelectricity to the total flexocoupling 
tensor as is clear from Eq. 0. requires the values of the elastic constants Cmiij which 
we also evaluated using the dynamical matrix, specihcally for perovskite cubic structure 
using the relationship 


kl 


1 


/ j ^ip^kp' 


( 12 ) 


p,p' 
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The technical details of the calculations are as follows. The calculations were done for 
cubic STO exploiting the Quantum ESPRESSO (QE) [H] ab initio package within the GGA 
PBE exchange-correlation functional with ultrasoft pseudopotentials [TB] . We have used an 
automatically generated uniform 16x16x16 grid of k-points, the kinetic energy cutoff for 
wavefunctions was 80 Ry. The calculations of phonon spectrum are done using Density 
Functional Perturbation Theory (DEPT) [16] as implemented in the PHonon code. The 
energy threshold for self-consistency was chosen to be equal to 10“^*^ Ry with the help of 
convergence tests for gamma point phonons. 

To hnd the independent components of tensors Xij = X^ijy ^ij = fijki, and Cijki (y, 

Af, /ii, fi2, /44, Cii, Ci2, C44), we evaluated the wave-vector dependence of the dynamical 
matrix Aip^iipi{q) for the [100] and [110] directions of the reciprocal space. This provided us 
with the components of the and matrices needed to finalize the calculations 

using Eqs. 


-O- 

The results of the calculations are presented in Tables |T] and [TTj A remarkable observation, 
which forms the main message of this Letter, is a very strong renormalization of the flex- 
ocoupling tensor by dynamic flexoelectricity: it is seen that the corresponding components 
of the tensors fijki and can drastically differ. This implies that the traditional neglect 
of dynamic flexoelectricity in the dynamic electromechanical simulations incorporating flex¬ 
oelectricity is inadmissable. Notably, the dynamic flexoelectricity can play an essential role 
in the frequency dependence of the electromechanical response close to the frequencies of 
mechanical resonances of the sample. The point is that such response is expected to be 
sensitive to the dynamic flexoelectricity only above the relevant resonance frequency [6]. 

The results of our calculations of the dynamical matrix provide us with an alternative 
way of evaluating some components of the material tensors of STO directly from the phonon 
dispersion curves. This will provide a cross-check of the above results. In addition, the 
method presented below will offer a way of evaluating these tensors from the experimental 
phonon spectra of the material. We suggest htting the long-wavelength part of the low- 
energy spectrum of the crystal to that obtained from the continuum theory, incorporating 
dynamic flexoelectricity. 

The suggested approach uses the spectra of transverse acoustic (TA) and soft-mode trans¬ 
verse optic (TO) branches for the high symmetry [100], [110], and [111] directions of wavevec- 
tor. In this case, the phonons are not accompanied by any wave of the electric held and the 
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M (xlO-8^) 

6 ±0.5 

hi (V) 

1.11 ±0.05 

fl2 (V) 

-1.30 ±0.05 

/44 (V) 

-0.28 ±0.05 

fit (V) 

-2.1 ±0.1 

fit (V) 

-2.5 ±0.1 

fit (V) 

-1.5 ±0.1 

x/eo 

2400 


TABLE I. Calculated static (/), dynamic (M), and total (/''“*) flexocoupling coefficients for cubic 
SrTiOs (a = 3.886 A) under a pressure of 78 kBar using dynamical matrix, x/eo is the relative 
dielectric susceptibility. 


(10^'^) 

(a) Dyn. mat. 

(b) Ph. disp. 

(c) Exp. 

Cll 

2.77 ±0.05 

3.1 ±0.1 

3.16 

Cl2 

1.06 ±0.05 

0.9 ±0.1 

1.03 

C44 

0.92 ±0.05 

1.0 ±0.1 

1.22 

|M|^ (V) 

1.2 ±0.1 

1.1 ±0.4 


\flt\ (V) 

1.5 ±0.1 

1.5 ±0.2 

1.2-2.2 

\flt - flt\ (V) 

0.2 ±0.1 

< 0.5 

1.2-1.4 


TABLE 11. Material parameters for cubic SrTiOa (a = 3.886 A) under a pressure of 78 kBar 
obtained using the dynamical matrix and from analysis of phonon dispersion curves, c is the 
stiffness tensor. is defined in Eq. (20). p is the density {p = 5174^) calculated by the 

mass of atoms in the unit cell divided by the cell volume, (a) Calculated from dynamical matrix, 
(b) Obtained from the simulated phonon dispersion curves, (c) Obtained from the experimental 
phonon dispersion curves [TTHIS]. 


continuum theory equations, Eqs. ([^ and ([^ with Ej = 0, yield the following dispersion 
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equation for the frequency u of the transverse modes [6] : 






2 CeffQ' 
= 


Uq = 


PI 

X~' + ^efrg^ 


P 7 

where the coefficients /eg, Ceg, and are dehned according to the direction; 


( 13 ) 

( 14 ) 


| 110 ] 


C44, /eg 9eS 5^44 

( 15 ) 

Ceg = ^(Cll — C12) 


/eg = lifll - /12) 

( 16 ) 

5 'efr = lidii — 912 ) 


Ceg = |(Cll — C12 + C44) 


/eg = |(/ll “ /l2 + / 44 ) > 

( 17 ) 

fi'eg = |(fi'll — 9l2 + fi'44) 



|111] 


and 7 can be found from the frequency of the lowest TO mode at T-point and the suscepti¬ 
bility, namely [HI [20l |2l] 

7 = (18) 

wo(0)x 

The roots of Eq. (13) with coefficients dehned by Eqs. (15) and dl?]) yield the frequencies 


of the corresponding double degenerate TA and TO modes. The roots of Eq. (13) with 


coefficients dehned by Eqs. (16) yields the frequencies of the TA and TO modes that are 


polarized along the [110] direction. 


One can perform the analysis of the dispersion of the acoustic branch by using Eq. (13) 


and get information on the parameters /eg and M. In the case of weak interaction between 
the branches (lowest order in g), the relative shift of the acoustic branch can be found as 


2 


2 2 
U — 




eg 


P7(a;2 - u\) ’ 


where is dehned as 


fX = /=« - . 

P 


( 19 ) 


( 20 ) 


consistent with Eq. (§. As is clear from Eq. ( [I^ , the deviation of the TA branch from 
the linear dispersion corresponds to the diherence — uj\ proportional to The 
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repulsive character of the interaction between the TA and TO branches (or lowering of the 


acoustic branch) is seen from the sign of rhs in expression (19). As for its magnitude, it is 


mainly controlled by the total flexocoupling coefficient (20), which has both dynamic and 
static contributions. One should note that, attempting to evaluate the flexocoupling tensor, 
such £t was done in several papers HEM], however neglecting dynamic flexoelectricity, 
which is unacceptable in view of the results presented above. Thus, the lowest order in the 
analysis of |a;^ —a;^| can only provide information about the total flexocoupling coefficient, 
not distinguishing between the static and dynamic contributions. The distinction between 
/ and can be obtained by expanding the solution to Eq. (|^ up to the next power in q: 


u 


= 


gVeff* - 




P7(w, 


OJ, 


( 21 ) 


A form of (21) suitable for analysis reads 






/eff - q'M 


x(/, 


eff 


( 22 ) 


It is clear that by fitting T(g^) to a linear function of one can get information on M 
and /gg*. However, the resulting information on M and the components of the flexocoupling 
tensor is limited. Only the absolute values of M, /n —/ 12 , and /44 can be provided from such 
analysis. We have checked that the analysis of the phonon spectrum for other directions 
than those addressed above unfortunately does not provide any additional information on 
the components of the flexocoupling tensor. Actually, all the information can be obtained 
from the dispersion curves for any two of the [100], [110], and [111] wavevector directions, 
while the third direction may serve as a crosscheck and an evaluation of the inaccuracy of 
the method. 

We have applied the described above method to evaluate the flexocoupling and flexo- 
dynamic tensors directly from the simulated phonon spectra. The phonon spectra were 
obtained for the [100], [110], and [111] wavevector directions on 10x10x10 automatic 
Monkhorst-Pack grid of q-points, in the framework already discussed above. The calcu¬ 
lations were performed under stabilizing pressures. An example of the TA and soft mode 
TO spectra are shown in Fig. [Tj The TA-TO repulsion effect caused by the flexoelectric 
coupling is clearly seen in this figure. To get information on the Mij and fijki tensors, the 












long-wavelength part of the spectra was calculated under a pressure of 78 kBar where the 
inter-mode coupling is a maximum. 



FIG. 1. Calculated phonon dispersion of transverse acoustic (TA) and soft-mode transverse optic 
(TO) branches with wavevector q = ^(C)0,0) for cubic SrTiOs stabilized by pressure. The red 
lines correspond to lower pressure {p = 78 kBar, a = 3.886 A) when the flexoelectric interaction 
between the TO and TA branches is strong. Here one observes lowering of the TA caused by 
flexoelectric coupling. The blue lines correspond to high pressure {p = 135 kBar, a = 3.85 A) 
when the interaction is weak in view of the large distance between the TA and TO branches. 
One can also observe lowering of the TO branch and change of dielectric susceptibility: for low 
pressure (red line), we have a;o(0) = 31cm“^ and x/^o = 2400, whereas for high pressure (blue 
line), a;o(0) = 107cm“^, and x/^o = 220. 


As a hrst step, M and were obtained from a linear £t for \I/(q'^) as a function of 
for the [100] and [110] directions. An example of such a fit is shown in Fig. It is 
worth mentioning that the reliable evaluation of the flexoelectric parameters directly from 
the spectrum is a challenging computational task. The problem is that, on one hand, one 
is interested in the long-wavelength limit of the spectrum, on the other hand, close to the 
F-point the violation of the acoustic sum rules corrupts the simulated spectrum (TU]. Thus, 
for this evaluation, only a limited interval of the eigenvectors should be used for the £t, 
as it was done in Figj^ Here, we used the simulated TA and TO frequencies for and 

Uq, respectively, p = 5174^. u\ and 7 were found using Eqs. (14) and (18) with the 
susceptibility and elastic constants obtained from the results of the same hrst principles 


calculations. Then, based on the results of the £t, using (20) we evaluated the values of 
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1 ^ 1 ) I/ll ~ /12I, and I/44I, which are given in Table [H} The spectrum for the [111] direction 
was found in a reasonable agreement with results obtained for those for the [100] and [110] 
directions, enabling us to evaluate the inaccuracy of our calculations given in this table. 





FIG. 2. Module of (22) (dots) and fitting line, which was used to obtain and M. To 

evaluate the frequencies of the transverse acoustic (TA) and soft-mode transverse optic (TO) 

branches entering in Eq. (22) were calculated for the [100] direction of wave vector, q = ^((",0,0), 
for cubic SrTiOs stabilized by pressure. 


Table [IT] gives a comparison between the results obtained with the two methods and with 
the available experimental data. Note that for the moment, experimental information is 
only available on some components of the tensor is available. It is seen from Table 
that the theoretical results yielded by the two methods are in good agreement with each 
other. The agreement between the theory and the data obtained from the Brillouin and 
neutron scattering experiments [MQ] is reasonable. 

The proposed method of evaluating the components of flexodynamic and flexocoupling 
tensors from the simulated spectra can also be applied for the determination of these compo¬ 
nents from the experimental phonon dispersion curves of materials exhibiting a ferroelectric 
soft mode. Thus, having information about phonon dispersion, dielectric susceptibility, and 
density of the material, one can exploit the above described method to hnd the flexody- 
manic tensor and certain combinations of flexocoupling coefficients experimentally. It is of 
importance to mention that, formally, this method can be applied to the phonon spectrum 
of any material, however in the absence of low-energy TO mode, it will not give reliable 
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information on the songht tensors. The point is that in such a situation the self-dispersion 
of the TA mode (which is always present and related to other effects) is expected to be of 
the same order of magnitude as the dispersion caused by flexoelectricity. 

It should be noted that for the moment the proposed method seems to be the only one 
enabling the evaluation of the strength of the static bulk flexoelectric effect (i.e. flexocoupling 
tensor) since static experiments in a hnite sample (like those done by Zubko et al [22]) yield 
the sum of the static bulk and surface contributions, which, in general, are comparable 

In conclusion, using ab initio DFPT calculations, for strontium titanate, we have evalu¬ 
ated the strength of the dynamic flexoelectric effect and compared it to that of the static bulk 
flexoelectric effect. We have found these effects to be of a comparable magnitude, in agree¬ 
ment with earlier order-of-magnitude estimates [12]. We identihed the situation where the 
dynamic effect is a few times stronger. Taking into account that, currently, all experimental 
studies and simulations involving flexoelectricity ignore the dynamic effect, we believe that 
our hndings are of importance. We have also offered a method for extracting information on 
flexodynamic and flexocoupling tensors using a simulated phonon spectrum. This method 
can be applied to the treatment of experimental phonon spectra, being currently the only 
method enabling an experimental evaluation of the size of the bulk flexoelectric effect. 

This project was supported by the grant of Swiss National Science Foundation, No. 
200020-144463/1, and by the grant of the government of the Russian Federation, No. 2012- 
220-03-434. 
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